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, We describe and evaluate a numerical solution strategy for simulating surface acoustic waves 

' (SAWs) through semiconductor devices with complex geometries. This multi-physics problem is of 

' particular relevance to the design of SAW-based quantum electronic devices. The mathematical 

^S) . model consists of two coupled partial differential equations for the elastic wave propagation and 

the electric field, respectively, in anisotropic piezoelectric media. These equations are discretized by 
the finite element method in space and by a finite difference method in time. The latter method 
yields a convenient numerical decoupling of the governing equations. We describe how a computer 
implementation can utilize the decoupling and, via object-oriented programming techniques reuse 
' independent codes for the Poisson equation and the linear time-dependent elasticity equation. First 

we apply the simulator to a simplified model problem for verifying the implementation, and there- 
after we show that the methodology is capable of simulating a real- world case from nanotechnology, 
involving SAWs in a geometrically non-trivial device made of Gallium Arsenide. 
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* S ' process of designing quantum electronic devices based on surface acoustic waves (SAWs) traversing piezo- 

electric media, it is necessary to determine the effect, on these waves of obstacles such as electrical gates on the 
r*^ ■ surface, and also the effect of the SAW on the low-dimensional quantum mechanical systems such as quantum wires 
and quantum dots. In general, the gates have a non-trivial geometry, which necessitate numerical simulation tools. 
The finite element method is well suited to handle complex geometries and is widely used to model piezoelectric 
devices Q, |^ H, |^ IM IB| - The effect on the low-dimensional quantum mechanical systems would be analyzed through 
J> , coupling the SAW simulator to both stationary and time-dependent Schrodinger equations 7j. This would require 
IT) ' the development of a fast SAW simulator but also flexible and portable code. 

C ! SAWs are modes of propagation of energy along the free surface of a material such that there is no decay along 
■ the direction of propagation but there is exponential decay into the bulk. SAWs have been used in electrical sensor 
technologies for many decades and have also been a useful tool in probing quantum electronic structures, for example, 
jy-^ quantum Hall liquids Q . Recently, much experimental work has been done in the field of acoustic charge transport 
whereby a SAW across a GaAs/AlGaAs heterostructure is used to capture a single electron and then transport it 
along a one-dimensional quantum wire This would be useful in developing an accurate current standard, but 
^ ' more challenging proposals to use this in the burgeoning field of quantum information processing have been proposed 
' [l^lT^ IT^ITslll^ . SAWs have also been utilized for both static quantum dot 15] and photo-luminescence experiments 
|l6lll7| . The time and resources required to build such devices are immense, and therefore the mathematical modelling 
I of these devices before the physical construction is advantageous. This approach requires the solution of the continuum 
Oh' electromechanical equations of motion in a piezoelectric medium. The method of partial waves can be used to 
obtain simple analytical expressions for the waves in the bulk material, but the the solutions say nothing about the 
effect of g ates on the surface. Attempts to solve the governing equations analytically for devices which do have surface 
gates |l9l l20l l2l| involve simplifications, and the accuracy of these approximations remains uncertain. 

There is a vast amount of literature on the three-dimensional finite element analysis of piezoelectric devices 



- - ' Q Hi Hi ^ IE 111 1 and these methods are exploited in the field of ultrasonics to design control systems involving 
piezoelectric actuators and sensors Commercial software such as ANSYS and ABAQUS can be used to simulate 
electromechanical phenomena but lack the fiexibility to couple to quantum mechanical calculations such as iterative 
Poisson-Schrodinger 0. Therefore it is often necessary to develop ones own computer code for modelling such devices. 

In this paper, we formulate and evaluate a finite element based solution method for the equations governing SAWs 
in piezoelectric media, and we formulate it in a flexible and portable manner to allow the ability to interface with 
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other code. The implementation is performed in an object-oriented style in order to incorporate existing solvers and 
to enhance portability of the code. Such a tool can be valuable in the design of micro- and nano-scale devices. We 
describe a set of boundary conditions that are capable of efficiently exciting SAWs and demonstrate the propagation 
of SAWs through a GaAs-based device. Although our applications are specific to GaAs, the formulation is general 
enough to allow simulations, with the same code (or with small modifications), of any crystal structure provided the 
elastic and piezoelectric material parameters are known. 

This paper is organized as follows. In Section ^ we give an overview of the mathematical model underlying 
piezoelectricity and precisely state the coupled partial differential equations we aim to solve. Sections IIIII and IIVI 
concern finite element formulations of these equations and an overview of the computational algorithm is described. 
In section [3 we describe the class based simulator approach used in our code design which simplifies implementation 
and increases reliability. Section IVII presents a verification of the implementation and an empirical estimation of 
convergence properties of the numerical solution method. A real-world application of the methodology, concerning a 
SAW problem in nanotechnology is the subject of Section IVlII before we make some concluding remarks in the final 
section. 



II. MODELLING OF PIEZOELECTRIC MATERIALS 



For piezoelectric materials, the electric displacement depends on both the applied electric field and mechanical strain, 
and the stresses depend on both the applied mechanical strain and applied electric field. The electric displacement 
Di is given by 

D, = efjEj+eijkejk, (1) 

where Ej denotes the electric field, is the permittivity tensor, e^fe is the piezoelectric coupling constant, £ij is the 
strain tensor. The stress tensor Oij is given by 

o-y = -SkijEk + cfjkl^kh (2) 

where Cijki represents the 4th rank tensor of elastic parameters. Summation over repeated indices is implied in this 
section, and the superscripts S and E denote that the quantities were measured under constant strain and constant 
electric field, respectively. 

In the absence of free charges within the material. Gauss's law requires the divergence of the electric displacement 
to vanish: 

V D^ Pfrcc - . (3) 

Combining Eqs. (Q) and |(2J), and the relation 

between the electric potential and the electric field together with the strain-displacement relation for small 
strains. 



where Ui is the mechanical displacement field, we can derive a scalar partial differential equation for 0: 



dxi ' dxk dxi dxk 

This equation couples the potential (j) and the displacements Ui in the medium, and can be used to compute the 
electric field. Equation I^J assumes the quasi-static approximation, which requires that very little energy is carried 
away by electromagnetic waves. 
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The displacement field in the medium is governed by Newton's 2nd law of motion combined with the appropriate 
constitutive law. The former equation reads 

S^^ = ^ + (7) 

where bi denotes body forces, and the double dot in Ui denotes a second-order partial derivative in time. As Eq. Q 
stands, it contains no damping term, which is irrelevant in our application setting. The constitutive law stated in 
Eq. (PI relates the stress tensor atj to the strain tensor e^, and Eq. ||SJ) relates £ij to the displacement field Ui. 
Combining Eqs. (0), 0, 01 we arrive at an equation for Ui in terms of 0: 



U uui U U(p 

= Q^c^^MQ^^ + 9b. + g^e.,,— . (8) 

From Eq. © we see how the mechanical motion is affected the electric field. In most materials, this coupling is weak, 
because Cijk ^ '^fjki^ the effect of the piezoelectric potential on the mechanical deformation is negligible. On 
the other hand, the effect of the mechanical deformation on the potential is always significant according to Eq. 0. 
However, in cases where an external electric field is applied through surface gates, the term in (j) may be sufficiently 
large to contribute significantly to the mechanical motion, and these cases receive the focus of the present paper. 

To illustrate the nature of these equations we explicitly write out the coupling terms for the piezoelectric material 
GaAs. The elastic parameters of a crystal are usually given in a coordinate system with its x, y and z axes parallel to 
the crystal X, Y and Z axes. For consistency with acoustoelectric charge transport experiments, we align the positive 
X axis with the crystal positive [Oil] axes, and the positive z axis with the crystal positive [100] axis. To achieve this 
we perform a 45 degree rotation of the crystal, around the z axis. 

In the transformed frame, the only non-zero components of the piezoelectric tensor are ei5 — —624 — 631 — —632. 
For generality, we illustrate our method using the different values of the piezoelectric tensor. Using the assumptions 
stated in the previous paragraph, and neglecting body forces, we can write out equations Eqs. © and © as 



d E dui d d(j) d d(j) 



d E dui ^ d d(j) ^ d d(j) 



(9) 



d E dui d d(j) d d(j) 

s d du, d du^ d du^, d du^ d duy d duy 

ox Ox Oy Oy Oz Ox Ox Oz Oz Oy Oy Oz 

In Refs. (2d) and l|2]T) analytical solutions of these equations were obtained, under the assumption that the cj) term 
in Eqs. l|^- Hll|) can be ignored. This means that the mechanical motion is decoupled from the electric field, but an 
electric field is induced from the mechanical motion. The reliability of such substantial simplifications is limited to 
cases where the external potential is small. Also, the obtained solutions are for two-dimensional cases only and 
are therefore of minor interest when studying the effect of surface gates. Solving the fully coupled system of PDEs 
demands numerical techniques like the one described in the next section. 



III. FINITE ELEMENT FORMULATION 



We shall use the finite element method in space and the finite difference method in time. The reasons for applying 
the finite element method are the need for handling geometrically complicated domains and the fact that the method 
works very well for elasticity problems without coupling to as well as for the Poisson equation for cj). 

We point out that the finite element formulation presented in this paper is independent of the choice of element 
shape, as this would be problem dependent. The elasticity part of the problem has three degrees of freedom - one for 
each component of the displacement vector. The electrical part has one degree of freedom for the electric potential. 
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A. The Poisson equation with piezoelectric coupling 

Before continuing further, it is convenient to introduce a superscript t to denote the time level, for example, (j)^ is at 
time level i. The electrostatic potential (j)^ due to mechanical displacements u\ in the piezoelectric material concerned 
is given by Eq. (|12|l . The finite element formulation of such a Poisson equation is well covered in the literature 
0,123) 13 113 ■ The basic idea is to approximate 4)^ by a Hnear combination of basis functions Ni, (f)^ k, cf)^ ~ ^j'Pj, 

insert 0^ in the Poisson equation, and demand the residual to be orthogonal to the space spanned by {A^i, . . . , Nn}. 
The Laplace term is integrated by parts. 

The only non-trivial aspect of the formulation in the present setting is the right-hand side, where the second-order 
derivatives of the elastic field demand integration by parts. The two first terms on the right-hand side of Eq. H12() 
give rise to an integral, which is straightforwardly integrated as 

, ^ . 9 dui d dui , ,^ 
ox ox ay ay 



dNi dul dNi dul , ,„ / du^ dui, 

'''^^ + ^^^^ 9^]^" " L ^^[^1^ 9^"- + e,,—ny]dT . (13) 

The terms in Eq. (|12|l . containing mixed derivatives, can be integrated by parts using a special form of Green's 
Theorem, 

di) f d(t> r 

dxdydz = / —ipdxdydz — <b (pipUzdV . (14) 



n dz ^ dz 



The terms give rise to the integrals 

f ,,,d dui d dui d dui d dui^^^ 
Jq dz dx dx dz dz dy dy dz 

. dN^dui dN^dui , dN.dul 9N,dul 

dz dx dx dz dz dy dy dz 

du^ du^ du'L du^ 

Ni[e3i—^nz + ei^—^Ux + 632^7— -|- e2i^—ny\dT . 
on ox dz dy dz 

The finite element method applied to the equation for (p transforms the PDE problem to a linear system of equations, 

= /, (16) 
where the 'stiffness' matrix P has its element given by 



P — 



dN, dNj dN, dNj dNi dNj 
dx dx dy dy dz dz 



dn . (17) 



The (p^ vector in Eq. (|16|) contains the values of (j)^ at the nodal points. The i-th component // of the right-hand 



side vector can be written as 



with 



<i efiV.^dr- f N,rdT+ / sdil, (18) 



dui dui dui dui dui dui 



y i_ y 

dx dz dy dz dx dy 



631^7— "2 + ei5^r-nx + e32^r^«z + 624^^^ + Cis^^Tl^ + 624 ^T""!;' 



dN.dui dNi dui dN.du^ dN.du^ dN, dui dN, dui 

s — 631— l-ei5— 1-632-^ \- £24-::: l-eis-::^ 7{ ^^24-;: . 

dz dx dx dz dz dy dy dz dx dx dy dy 

When the displacement field entering is known, any standard Poisson solver can be used to compute cj/ 



5 



B. The elasticity problem with electric field loading 



For the finite element formulation of Eq. ((Sjl it is easier to first start with Eq. Q and insert the constitutive law 
given by Eq. |(2Jl and the strain-displacement relation given by Eq. in the finite element integrals. The time 
derivative in Eq. © can be approximated by a second-order accurate finite difference. Sampling Eq. (jHJ at time level 
£ then yields 



^ aT^ -9^ + ^^- ^^^^ 

where At is the time step length and quantities with the £ superscript are functions of space only. This time 
discretization introduces an operator splitting such that the originally coupled governing equations can be solved in 
sequence. No accuracy is lost by this operator splitting beyond that implied by the finite difference itself in Eq. l(T^ . 
More precisely, afj contains (j)^, which is known, such that we can easily solve for m^"*"^ using old values of cj). On the 
next time level {i + 1), we can find 0^+^ using the recently computed wf"*"^. 

The main motivation for the explicit time differencing in Eq. (|19|1 is numerical efficiency: (i) we decouple the 
equations, which simplifies the numerics and the implementation, and (ii) there is no need to solve large sparse linear 
systems of equations in the elasticity part of the problem (if the mass matrix is lumped). Nevertheless, the time 
difference Eq. leads to a conditionally stable scheme, where At must be of the order of the smallest element size. 
In wave propagation problems, uniform high resolution is frequently needed in space and time, typically compatible 
with the stability restriction, which makes such explicit schemes appropriate. On the other hand, in applications 
where adaptive grids with great variation in element size are needed, one may benefit from implicit schemes (for 
example, of Newmark type) 24] . 

The finite element formulation of Eq. (|19|) is easiest to express if we switch to a typical "finite element engineering" 
notation, especially when we deal with anisotropic media. The stress and strain tensors are expressed as vectors. 



^ XX T ^yy T ^ zz T '^^yz : zx T '^^yx ) 



T 



The constitutive law can then be written as 



^De'+p', (20) 



where the p term represents the loading from the electric field, and Z) is a symmetric 6x6 matrix of the elasticity 



coefficients (previously denoted by c^^j): 



Cll C12 Ci3 Ci4 Ci5 C16 
C12 Ci3 Ci4 Ci5 C16 



D 



Cl3 Ci4 Ci5 C16 



Ci4 Ci5 C16 



(21) 



Cl5 C16 



C16 

The displacement field at a time level £ is approximated according to 



n 

where is the value of the displacement field at node j at time level I. The strain-displacement relation then becomes 



6 
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TV- 



































(22) 



The comma notation here denotes partial derivative: Ni_x = dNi/dx. 

A finite element formulation of Eq. p9|) using the aforementioned notation becomes (see Ref. !23 for a detailed 
derivation without the (j) term) 



n n n on 



(23) 



Here, is the traction on the boundary, arising from integrating the divergence of the stress in Eq. H19(l by parts. 
Inserting the constitutive law in the term J crd^l yields 



^ j j BjDBjdn j + J Bfp'^dil . 



The final discrete equations arising from the equation of motion can be written as 



u 



-1 = 2u' - u 



£-1 



AfM i-Ku'^ +(3'^ + 



(24) 



where M is an efficient inverse of the mass matrix M resulting from the gNiNjdfl integral. We construct M 
as the inverse of the lumped mass matrix (using the row-sum technique to lump the matrix). The matrix K stems 
from the standard "stiffness" term BJ DBjdfl representing anisotropic elasticity, is the effect of body forces 
and surface tractions, and is the contribution from the electric field. The i-th block (arising from node i) in 
takes the form $j — B^ p^, which from Eqs. ^ to ((TT|l results in 



-InN.[fye.2^ + £e,,^]dn 



dy dz 



dNj dip 



dz dy JaO^'U'-J^ Oz '"y ' "--^4 Qy 



/o[ei5l^^ + + Ion mei,^n^ + e,,^ny]dr 



(25) 



Our governing equations require initial conditions. Let us assume that the elastic body is at rest such that dui/dt = 
0. The external electric field is then turned on. After the initial transients the displacement field have faded out, and 
we have a stationary initial state of our system, modeled by the equations 



P4>" = /"(«"), 



(26) 
(27) 
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These two coupled equations are solved by an iterative Gauss-Seidel-like technique, i.e., the equations are solved one 
at a time, using the most recent approximation of the other field in the right-hand side term. Each linear system is 
solved by a MILU preconditioned conjugate gradient method [2^. The solution of Eqs. and (P7|l . along with the 
assumption of stationarity, du/dt — 0, comprise the initial condition. 

As boundary conditions, we either have prescribed traction components or prescribed displacement components, 
along with prescribed electric potential or prescribed normal component of the electric field. The displacement 
equation needs three boundary conditions at each point at the boundary, while the equation for needs one condition 
at each point. 

From du/dt = at t = it follows by a second-order difference approximation that = u^^. From Eqs. H24|l and 
(|27|l it follows that — — . The loads removed at t = 0+ will first come into play at the second time level. 
The computational algorithm can now be formulated as follows: 

i — Q (time level counter) 
A; = (iteration counter) 
while e < Ecrit 

solve Eq. ^ w.r.t 

solve Eq. l|77|l w.r.t u^'^ 

compute e =|1 u°^^ - u°^^~^ \\ + || cfP'^ - cj)"'''-^ \\ 
fc «— fc + 1 end while 
^1 = = m" 

for i — 1,2,3,... until end of simulation 
solve Eq. w.r.t 
solve Eq. (Uni w.r.t 0^+^ 

Comments on Stability. The stability criterion of the explicit finite difference scheme in time, when decoupled from 
the electric field problem, requires At < 2/uJma.x 27], where o^max is the highest natural frequency of the vibrating 
system. One may find Wmax as the square root of the largest eigenvalue of the problem K — XM — 0. An approximate 
bound on Wmax can be estimated from a relation Wmax = 2c//ioff: where c is the speed of elastic waves and /i^ff is the 
smallest effective element length. The stability criterion now reduces to the common Courant, Friedrichs, and Lewy 
(CFL) condition: 

M<aKs/c. (28) 

Here, a is a factor to be adjusted since the hcs parameter is normally roughly computed from element sizes or 
application of Gerschgorin's theorem applied to the underlying eigenvalue problem. For wave problems the CFL 
condition is frequently not particularly restrictive since the length and period of a wave are usually proportional, 
leading to a natural choice oi At/h — const. 



IV. AN OPTIMISED FINITE ELEMENT FORMULATION 



We can speed up the numerical computation by replacing the complete assembly of the right-hand side vectors in 
Eqs. (|15|l and (|24|l by a matrix-vector product. For example, consider the contribution to the vector / in Eq. Hlt)|) 
due to the coupling from the mechanical motion (omitting the superscript £ for the remainder of this section): 

f dNi dux dNi dux dNi duy dNi duy dNi du^ dNi du^ ,^ 

/ ^31^^ + '=15^^ 7^ + e32^r 77^ + 624^; 77^ + 615^^ ^+624^^ H"^^- ^9) 

Jf2 oz OX OX az az oy ay az ax ox ay ay 

Inserting the finite element expansion Ux ~ X]j=i ^j^^-^ii ^"^^ similar expansions for Uy and u^, the expression in 
Eq. (123) results in 



f ON, ON, (x) dN,dN. 



dx dz ^ 



ez2-^ 7\ — 

oz Oy ■' 



(30) 



Oy Oz ■' 



ei5^- 5— w,- 

ox Ox 



624- 



dy dy 



■uT'dn. 
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This expression can be written as a matrix vector product Lm, where the matrix L consists of n x n blocks, each 
block of size 1x3. For the coupling of node i and j, the 3x1 block looks like 

dNidNj dNidNj dN.dNj dN, dNj dN.dNj dN.dNj 

631-?^ T\ 1-615-;::; 7: — '632-^ ^ 624-?; ^ — '^15-- [-624-?; ^ — , (31) 

oz OX OX Oz Oz Oy oy Oz Ox Ox Oy Oy 

The Poisson equation - Eq. for can now be written as 

P<^^+i = Lw^+i . (32) 

At each time level we can hence avoid the costly finite element assembly process, since P and L are constant in time 
and the right-hand side of the linear system is obtained by an eflticient matrix- vector product. This may result in a 
significant speed-up of the solver. For large number of unknowns, the speed up is experienced only if we use a method 
of complexity or order n, like multigrid, to solve Eq. (|32|l . 

Equation 124(1 is already on a favorable matrix- vector algebra form, except for the coupling term Examining 
Eq. H25() . we realize that this matrix can be written as — Q4>^, where Q is a time- independent n x n matrix of 
3x1 blocks. The contribution to block i in is made of the coupling Q^j between node i and j and the block j in 
4>^: *f = Qy 0^^, where 



Q^ 



r dNidNj , dNidNj r r dNj . dNj i 

r dNidNj . dNi dNj . r j,jr dNj , dNj i 



Now the ^•'^ term in Eq. I|24|l can be computed as a matrix-vector product Qcj)^, and no assembly process is required 
to construct Eq. H24|l at each time level. In our algorithm we must run assembly processes at i — to construct the 

matrices M , K, Q, L, and P. 



V. OBJECT-ORIENTED IMPLEMENTATION 



A significant trend in modern software development is to formulate numerical algorithms such that reliable and 
well-tested software components can be combined together to form a new simulator. Our numerical approach was in 
particular inspired by such an approach. With the time discretization we were able to split the coupled Ui-(j> system 
such that at each time level we first solve for a new displacement field (ui) and then we solve for the corresponding 
electric potential (0). In each of the two equations, the effect of the other only enters through a right-hand side 
"forcing" term. This allows us, in principle, to reuse a solver for elastic vibrations and a Poisson solver, as long as 
these solvers can implement a new right-hand side term that couples to another solver. From a principal point of 
view, this idea is simple and attractive. However, the implementation of the idea in practice may be less feasible if 
the design of the underlying solvers is not sufficiently fiexible. 

To allow for building the compound Ui-(f) solver from separate Ui and (j) solvers we propose to apply principles 
from object-oriented programming. The time-dependent elasticity solver and the Poisson solver are realized as two 
independent objects, implemented via classes in C-I--I-. Actually, we have reused the Poisson2 class from Ref. 23 as 
the Poisson solver, and the ElasticVibl class from Ref. l23l as the elastic vibration solver. The latter is a subclass of 
Elasticity2, a pure quasi time-dependent elasticity solver (without the acceleration term in the momentum equation). 
These classes are implemented using building blocks from the Diffpack library . 

The Poisson2 and ElasticVibl classes have no knowledge of each other. The only common feature is their design. 
This design is crucial for reuse of the classes to solve the coupled problem, but the design approach suggested for 
Diffpack solvers [23112^12^ has proven to be successful in this respect. Diffpack solvers are realized as classes containing 
objects for grids, scalar/vector fields, linear systems, etc. The coefficients in a PDF are evaluated through virtual 
functions. When the solver is stand-alone, these virtual functions contain mathematical expressions or measured data, 
but when solvers are combined, the virtual functions are reimplemented in subclasses and connected to data in other 
solvers. This principle reffects the underlying mathematics: the coefficients in a single PDF are considered known, 
but in a system of PDEs, the coefficients typically couple to unknown quantities governed by the PDEs. The present 
coupled system is an example where the right-hand side in the Poisson equation couples to the primary unknown 
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in the elastic vibration solver, while a right-hand side term in the elastic vibration solver couples to the primary 
unknown in the Poisson equation. When the solvers are used independently there are no such couplings. 

Figure n] shows an outline of the class design for the compound solver. From class Poisson2 we derive a subclass 
Poisson2_glue, which reimplements the virtual function in Poisson2 for evaluating the right-hand side in that equation. 
In this new function we need to compute an expression involving the displacement field available in the ElasticVibl 
class. A manager class, simply called Manager, holds pointers to all classes for the individual PDEs in the system 
of PDEs, and each PDE class holds a pointer to the manager class, thus enabling a two-way communication. With 
these pointers, we can connect data from any other solver class to the Poisson2_glue class. In the virtual function 
evaluating the right-hand side we can typically call mng->elastic->u to get a vector field object for Ui that we can 
evaluate at the current integration point. This object has support for evaluating derivatives of the field as well. 

Similarly, we derive a subclass ElasticVibl_glue of ElasticVibl, add a pointer to the manager class, and override 
the virtual function for evaluating the right-hand side. Now we need to compute expressions involving cf), but this is 
easily accomplished by the pointers, e.g., mng->poisson->u if u is the name of the unknown scalar field in the Poisson2 
solver. 

The "glue" classes Poisson2_glue and ElasticVibl_glue are very small compared to the real solver classes they inherit 
from. A "glue" class typically needs about a page of code unless it adds additional computations. The manager class 
is a bit more comprehensive since it implements the overall solution algorithm, i.e., the time stepping and the calls 
to the independent solvers. 

The benefit from using object-oriented programming in the way we have outlined is that the original well-tested 
Poisson and elastic vibration solvers can be reused in a system of PDEs without any modifications. The original solver 
classes reside in separate files and are hence not subject to any side effects from editing the source code. The "glue" 
class is also in a separate file and enables the original solver to speak to a manager in charge of solving a compound 
system of PDEs. The advantages are clear: reliability is increased by reusing well-tested solvers, and the compound 
solver is modularized. Our experience is that this design reduces the development time significantly. Especially the 
debugging phase is greatly simplified. At any time, the underlying solvers can be trivially pulled out of the compound 
system and verified independently. 



Elasticity2 



Poisson2 



ElastVibl 



Poisson2_glue 



ElastVibl_gliie 



poisson elastic / 

mng \ / ' mng 



Manager 



♦ is-a relationship 

* has-a relationship 



FIG. 1: Relationships between classes. The base classes, Poisson2 and Elasticity2 perform the solution of the standard 
time-independent Poisson and elasticity equations, respectively. ElasticVibl is derived from Elastcity2 and solves for time- 
dependent elastic motion. The classes Poisson2_glue and ElasticVibl.glue implement the coupling between the electrostatic 
and mechanical equations, and these classes are controlled by the Manager class. 
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VI. VERIFICATION 



Before showing numerical results for a physically relevant application, we report on the estimated numerical accuracy 
of the code, thereby establishing evidence for the correctness of the implementation. To investigate numerical errors, 
we compare numerical results with an exact solution. The exact solution is based on the assumption of a displacement 
field in z direction only, depending on x and t only. Physically, a normal traction is applied at all x — const and 
y — const boundaries to avoid displacements in the x and y directions, but in the code we implement this situation 
by essential boundary conditions Ux = Uy = 0. We also assume that (j) = ^^'^ that the material is isotropic. 

The governing equations then reduce to 



These equations can be scaled to yield 



= ^^^+"7^2' (35) 



dx"^ dx^ 
dx^ dx"^ 



(36) 



The constant a is zero or unity corresponding to whether the elasticity problem couples to the electric field problem 
or not. (The three coefficients in the original system are scaled away by choosing an appropriate time scale, (j> scale, 
and Ui scale.) One possible solution of Eqs. (|35|1 - H3()|) reads 



Uz{x, t) — (f>{x, t) — cos{k\/l + at) sin(A:a::) . (37) 

In our tests we choose $7 as a three-dimensional beam and k = with n being an integer and L the length of the 
beam. The ends of the beam are then fixed. 

Our verification procedure consists of estimating the convergence rate of the numerical method. A scalar error 
measure e is expected to behave like 

e = Ah'' + BAt\ 

where A and B are constants independent of the element size h and the time step At. From the involved approxi- 
mations, we expect s = 2 and r = 1 + q, where q is the order of the polynomials in an element. In particular, for 
linear or trilinear elements, r = 2, and if h and At are chosen such that h/At = const, we see that e ^ . From two 
successive experiments, (/ii_i,ei_i) and (hi,ei) we may estimate a convergence rate from 

lnej/ei_i 



In hi /hi 



We have investigated three error measures, given as different norms of the error field E ^ Uz — Uz, where iiz is the 
numerical solution and Uz is the exact solution: 



miL^n) = I \E\dn, \\E\\L2(n) = ( / E^dnf^, \\E\\L^in) = sup(|i?(x)|, x e f!) . 

Table I displays the values of the error norms and the associated estimates of the convergence rates r^. As the grid 
spacing h and the time step At are reduced, the numerical solution converges to the exact solution with the expected 
rate of two. This provides evidence for the correctness of the implementation and indicates that our overall solution 
method, including the operator splitting, is of second order in time and space. The spatial convergence rate is expected 
to increase with the order of the polynomials used in the elements. 



VII. A SURFACE ACOUSTIC WAVE APPLICATION 



In order to show that our suggested numerical model may be applied to real-world phenomena, we apply it to a 
problem in quantum electronics, specifically in acoustic charge transport. We simulate a SAW propagating through 
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TABLE L Error norms and estimated convergence rates for waves in a piezoelectric material, with h — lOAt. 

a nano-scale substrate of GaAs. SAWs are particular solutions of Eqs. © - H12|l such that they propagate without 
decay on the surface of a material but decay exponentially into the bulk. The solutions typically have the form, 

u, = [/,_j-e-'^«^-^e*('==^-'^*\ (38) 
i 

j 

assuming x is the direction of propagation, z is the direction into the bulk (z — > — oo), and the decay constants q, , 
may be complex, allowing for oscillatory decay (in an exponential envelope) into the bulk as is the case for GaAs [33 ■ 
The Ui^j and $j are constant amplitudes, and fe is a wavenumber. It is a straightforward mathematical procedure to 
determine the decay constants and wave velocity [T^ . However, we are interested in the more complicated dynamics 
taking place when an obstacle in the form of a charged metallic gate is placed on the surface. 

In acoustic charge transport, the time-varying electric potential accompanying the SAW is used to transport elec- 
trons which are trapped in the minima of the waves. To perform further manipulation on these electrons, for example, 
to remove some from a minimum, external electric fields must be applied, and this is achieved by applying voltages to 
metallic gates placed on the surface. In general, the gates have different electrical and mechanical properties to the 
bulk piezoelectric material and so we may expect to see interesting effects if for example a SAW is passed through the 
compound structure. Here, we perform three simulations where we pass a SAW through a piece of GaAs. The first 
will be a bare SAW, as a check to see we do excite the required modes. The following two will have a metallic gate of 
contrasting mechanical properties to GaAs with a static voltage applied on it. In the first of these, we decouple the 
mechanical motion from the electric field (while keeping the electric field coupled to the mechanical motion). In the 
second we allow the full mutual coupling between the electrical and mechanical fields to take place. The dimensions 
of the gate are 600 nm x 400 nm x 200 nm. To magnify the effect of the compound material structure we use a 
fictitious material for the gate, similar in crystal structure to GaAs, but with the elastic constants and mass density 
changed by one order of magnitude. Moreover, we apply a relatively large voltage of 1.5 V so that the coupling of the 
mechanical displacement to the electric field is evident. 

In the numerical simulations we have implemented 8-noded hexahedral brick elements with 2x2x2 Gauss-Legendre 
integration. The spacing along the direction of the SAW propagation - Ax is chosen to be 100 nm. The time integration 
parameter At is 1.25 x 10"^ ns. 

A. Experiments with surface acoustic waves 

In realistic SAW-based devices, bulk waves and SAWs are excited by the application of a microwave signal to an 
interdigitated transducer. The waves propagate outward from the source, with bulk waves dissipating energy into the 
bulk and surface acoustic waves travelling through the free surface. In acoustoelectric charge transport experiments, 
the quantum devices are located several thousand microns away from the source where one would expect to see 
only surface acoustic waves. To simulate the entire system would therefore require vast computational resources 
(the wavelength of a SAW is one micron and approximately ten nodes would be required across one wavelength for 
a satisfactory description). Instead we describe a system of boundary conditions which excite acoustic waves with 
SAWs dominating over the bulk waves within a few microns of the source. 
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1. Excitation of SAWs. 

We point out that to be consistent with the formulation presented in the previous sections and hence acoustoelectric 
charge transport experiments, in our simulations the SAW travels along the crystal positive [Oil] axis, with the positive 
z axis aligned with the crystal positive [100] axis. 




FIG. 2: A schematic diagram of the experimental setup for simulating SAWs. The gate is centered at (8.9- 10 nm, 2.6-10^ nm) 

In order to excite the SAW modes, we apply a time-dependent Dirichlet boundary condition to the z component of 
the displacement field on a small region on the grid. This can be expressed as 

Ux — Uy — 0, Uz — Asm{2Tr ft), (39) 
xq < X < xi, z — 8000 nm . 

The thickness of the membrane is defined as xi — xq — 100 nm. To reduce diffraction effects from the surfaces - 
y = nm and y = 5200 nm we have forced Mj, = on those faces. The frequency / used is 2.7 GHz and the amplitude 
A is chosen so that the SAW amplitude is approximately 20 mV at a depth of 100 nm. Traction free boundary 
conditions for the mechanical equations and zero normal electrical displacements are implemented at the surfaces. 
We found this system of boundary conditions to be a particularly efficient means of generating coherent SAWs. The 
setup is shown in Fig. |21 

Figure |3] is a result of a simulation showing the ^ phase difference between the x and z components of the dis- 
placement vector, and the larger z amplitude, suggesting elliptical polarization in the sagittal plane. Figures^a) and 
EJb) show two-dimensional slices, in perpendicular planes, of the full three-dimensional solutions. The first shows 
the non-decaying nature of these waves along the propagation direction. The SAW wavelength can be seen to be 
approximately 1000 nm and its velocity is computed to be approximately 2700 ±20 ms~^, which is in good agreement 
with analytical calculations for the wave velocity. The second figure shows the surface nature of these waves; the 
greatest amplitude is near the surface and there is no observable decay in the direction of propagation but they decay 
exponentially into the bulk. On the far left of each image (close to where the boundary condition is applied), bulk 
waves may be observed but as the waves propagate towards the right, they dissipate all their energy into the bulk and 
the only remaining waves are the surface waves. Figure |S1 shows the SAW amplitude as a function of depth. We see 
that the amplitude undergoes oscillatory (complex exponential) decay into the bulk and is negligible a few microns 
below the surface, again in good agreement with analytical expressions derivable from Eq. |30|. 

2. The effect of the compound mechanical structure. 

Figure Ela) shows the magnitude of the displacement field as the SAW moving from the left of the figure to the 
right, passes through the gate when the mechanical motion is decoupled from the electric field. The gate is centered 
at (8.9 • 10^ nm, 2.6 • 10'^ nm). We observe a peak to the left of the gate (due to a reflection) and a trough (due 
to damping of the wave) to the right. We also just barely see vibrations moving away from the gate along the y 
direction. Since in this simulation, the mechanical motion is independent of the electric fields, this damping is due 
purely to the presence of a mechanical structure on the surface. The damping of the mechanical amplitude may be 
large enough to affect the electric potential. In this simulation, the electric potential is not affected significantly as 
shown in Fig. |H|[b). The central blur is the electric potential due to the gate. The SAW peak emerging from the gate 
has been damped in the central region. 
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FIG. 3: The relative amplitudes and phases of the displacements and Uz, parallel and perpendicular respectively, to the 
SAW propagation. 



3. The effect of electromechanical coupling. 

Figure [7Ja) shows the magnitude of the displacement field, at time 2 ns, as the SAW passes through the gate, in 
the case where mutual coupling between the electric and mechanical fields is allowed. As in the decoupled case, we 
observe a peak to the left of the gate (due to reflection of the wave) and a trough (due to damping of the wave) to 
the right. Again, we also see vibrations moving away from the gate along the y direction. Comparing Figs.EIa) and 
El^a), we see that there is additional mechanical deformation throughout the material due to the presence of a charged 
metallic gate. Also, the scale of Fig.[2fa) is shifted up compared to that of Fig.E^a). By looking at the displacement 
field at time i.e. before the arrival of the SAW, we can see the displacement resulting from the applied electric field. 
This is shown in Fig. Efb). 
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FIG. 4: (a) Plane waves of the electric potential, 100 nm below the surface, at time 1.6 ns. (b) Surface nature of the electric 
potential, at time 1.6 ns. The SAWs have a significant amplitude between z — 8000 nm and z = 7000 nm and have a significantly 
reduced amplitude below z — 7000 nm 
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FIG. 5: A typical curve of the electric potential as a function of depth into the bulk. The precise shape of the curve clearly 
depends on the exact time and spatial coordinate the data was taken as the amplitude of the SAW varies from —20 mV to 
+20 mV. However, all such curves have the oscillatory decay into the bulk becoming negligible within a few microns. 
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FIG. 6: (a) The magnitude of the displacement field U — ^UrUr , at time 2 ns, as the SAW moving from left to right, passes 
through the gate centered at (8.9 • 10^ nm, 2.6 • 10^ nm) and shown by a white square. The mechanical motion is decoupled 
from the electric field, (b) The total electric potential, at time 2 ns, as the SAW passes through the gate when the mechanical 
motion is decoupled from the electric field. The scale has been restricted so that the SAW potential is visible. The central 
blur, which is due to the large external potential has the value of 1.5 V at its center. The white square shows the position of 
the gate. 
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FIG. 7: (a) The magnitude of the displacement field U — ^UrUr , at time 2 ns, as the SAW moving from left to right, passes 
through the gate centered at (8.9 • 10'' nm, 2.6 • 10^ nm) and shown by a white square, when full coupling between electric and 
mechanical fields is allowed, (b) The magnitude of the displacement field U = ^UrUr, at time 0, demonstrating the mechanical 
strains caused purely by the gate with a 1.5 V applied voltage centered at (8.9 • 10'' nm, 2.6 • 10^ nm) and shown by a white 
square, when full coupling between electric and mechanical fields are allowed. 
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VIII. SUMMARY AND CONCLUDING REMARKS 

Our aim with this paper is to suggest a computationally fast method for simulating SAWs in piezoelectric devices 

where stress, deformation, and a quasi-static electric field are fully coupled. The basic idea of the numerical scheme is 
to use a finite difference approximation in time that decouples the elasticity and the electric field problems such that 
these can be solved in sequence at each time level. This decoupling can also be explored in computer implementations, 
because independent solvers for anisotropic elasticity problems and Poisson problems can be joined together. We 
showed in particular how object-oriented programming techniques can realize such couplings in a very convenient 
way. 

Through numerical experiments in a physically one-dimensional wave propagation problem we have verified that 
the three-dimensional code reproduces the expected quadratic convergence in space and time if linear or trilinear 
elements are uscxl. Finally, we have applied the proposed numerical methodology to real SAW phenomena in a device 
with complicated geometry. We have described a system of boundary conditions which are capable of exciting SAW 
modes in a small computational domain. We have performed simulations where we have demonstrated the effects 
on the SAW of a charged metallic gate. The results indicate that the method is capable of predicting the expected 
complex elasto-electric dynamics in such a device. 

The methodology could easily be applied to simulate SAWs through devices such as GaAs/AlGaAs heterostructures 
with more complicated surface gate geometries. 

The principal limitation of the equation splitting is a stability criterion on the time step length. Implicit methods 
may remove time step restrictions, but at a cost of the need to solve large coupled linear systems at each time level. The 
computer implementation is also more involved and requires a special-purpose solution rather than just joining two 
well-tested equation components. However, for wave propagation problems one usually needs a fine mesh and a small 
time step to resolve the waves, typically leading to /i/At = const, which is in accordance with the stability criterion. 
Explicit time stepping approaches are therefore highly relevant and allow efficient algorithms and implementations. 

Forthcoming numerical work will focus on domain decomposition methods for parallelizing the solver and thereby 
enable simulation of large-scale SAW problems. The described time stepping approach and associated equation 
splitting are particularly well suited for parallel computing, because the elasticity problem can be made "perfectly 
parallel" , and very efficient parallelization strategies exist for the Poisson equation. 
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